% solver is a matlab function who do

% porosity upscaling - arithmetic mean;
% permeability upscaling - pythagorean means and flow based method:
% and create  the new IMEX .dat file with imex_export function. (if you
% want u can run it already to but now on test phase i disable this)

% by jose augusto f. ramos - last update date 22/08/11
% contact zeaugustofr@gmail.com

function solver(arquivo,outfolder,method,refino,transmi,i,j,k,ni,nj,nk,np,cells_i,cells_j,cells_k,...
            nDI,nDJ,nDK,PHI,KI,KJ,KK,DI,DJ,DK,wells)
        
% creating case file
outfolder = strcat(outfolder,num2str(ni),'_',num2str(nj),...
    '_',num2str(nk)); %an especific case
if refino == 1
    outfolder = strcat(outfolder,'\refino'); %refinament especific
end
if transmi == 1
    outfolder = strcat(outfolder,'\transmi'); %transmissibility especific
end
mkdir(outfolder); %create the output folder if didn't exist

% first - porosity upscaling
nPHI = porosity_upscaling(PHI,nj,cells_j,ni,cells_i,nk,cells_k,np);

% second - permeability upscaling
% now its depends on wich method was chosen. method is a vector with holds
% 4 numbers if one those is true the upscaling will happen. the program
% allows to make diferents upscaling at the same time.

% to keep the final .dat much cleaner we work with inlcudes that facilites
% the working a lot for the writing the new permeability e porosity matrix

nPHI_include = strcat(outfolder,'\nPHI.inc');

if method(1) == 1

    arithmetic_i = arithmetic_solver(KI,ni,cells_i,nj,cells_j,nk,cells_k,np);
    arithmetic_j = arithmetic_solver(KJ,ni,cells_i,nj,cells_j,nk,cells_k,np);
    arithmetic_k = arithmetic_solver(KK,ni,cells_i,nj,cells_j,nk,cells_k,np);
    
    arithmetic_output = strcat(outfolder,'\arithmetic_',num2str(ni),'_',...
        num2str(nj),'_',num2str(nk),'.dat');
    
    if refino == 0
        
        imex_export_3d(arquivo,outfolder,arithmetic_output,ni,nj,nk,nDI,nDJ,nDK,nPHI,...
            arithmetic_i,arithmetic_j,arithmetic_k,nPHI_include);
        
    else
        
        imex_export_3d_ref(arquivo,outfolder,arithmetic_output,i,j,k,ni,nj,nk,cells_i,cells_j,cells_k,...
            nDI,nDJ,nDK,nPHI,PHI,KI,KJ,KK,wells,...
            arithmetic_i,arithmetic_j,arithmetic_k,nPHI_include)

    end
    
    % %EXECUTE FILE IN IMEX (to do so discoment this)
    % aux = ('"');
    % aux1 = (' -wd '); 
    % aux2 = (' -f ');
    % system([imex aux1 outfolder aux2 aux arithmetic_output aux '&'])
    
end
if method(2) == 1
    
    geometric_i = geometric_solver(KI,ni,cells_i,nj,cells_j,nk,cells_k,np);
    geometric_j = geometric_solver(KJ,ni,cells_i,nj,cells_j,nk,cells_k,np);
    geometric_k = geometric_solver(KK,ni,cells_i,nj,cells_j,nk,cells_k,np);
    
    geometric_output = strcat(outfolder,'\geometric_',num2str(ni),'_',...
        num2str(nj),'_',num2str(nk),'.dat');
    
    if refino == 0
        
        imex_export_3d(arquivo,outfolder,geometric_output,ni,nj,nk,nDI,nDJ,nDK,nPHI,...
            geometric_i,geometric_j,geometric_k,nPHI_include);
        
    else
        
        imex_export_3d_ref(arquivo,outfolder,geometric_output,i,j,k,ni,nj,nk,cells_i,cells_j,cells_k,...
            nDI,nDJ,nDK,nPHI,PHI,KI,KJ,KK,...
            geometric_i,geometric_j,geometric_k,nPHI_include)
        
    end
    
    %     system([imex aux1 outfolder aux2 aux geometric_output aux '&'])
    
end    
if method(3) == 1
    
    harmonic_i = harmonic_solver(KI,ni,cells_i,nj,cells_j,nk,cells_k,np);
    harmonic_j = harmonic_solver(KJ,ni,cells_i,nj,cells_j,nk,cells_k,np);
    harmonic_k = harmonic_solver(KK,ni,cells_i,nj,cells_j,nk,cells_k,np);
    
    harmonic_output = strcat(outfolder,'\harmonic_',num2str(ni),'_',...
        num2str(nj),'_',num2str(nk),'.dat');
    
    if refino == 0
        
        imex_export_3d(arquivo,outfolder,harmonic_output,ni,nj,nk,nDI,nDJ,nDK,nPHI,...
            harmonic_i,harmonic_j,harmonic_k,nPHI_include);
        
    else
        
        imex_export_3d_ref(arquivo,outfolder,harmonic_output,i,j,k,ni,nj,nk,cells_i,cells_j,cells_k,...
            nDI,nDJ,nDK,nPHI,PHI,KI,KJ,KK,...
            harmonic_i,harmonic_j,harmonic_k,nPHI_include)
        
    end
    
    %     system([imex aux1 outfolder aux2 aux harmonic_output aux '&'])
    
end    
if method(4) == 1
    
    if transmi == 0 
        
    %permeability
    [fbasedk_i,fbasedk_j,fbasedk_k] = ...
        fbasedk_solver_3d(ni,cells_i,nj,cells_j,nk,cells_k,np,DI,DJ,DK,...
        KI,KJ,KK,transmi);

    else    
        
    %transmissibility
    
    TI = DJ*DK*KI/DI;
    TJ = DI*DK*KJ/DJ;
    TK = DI*DJ*KK/DK;
    
    [fbasedk_i,fbasedk_j,fbasedk_k] = ...
        fbasedk_solver_3d(ni,cells_i,nj,cells_j,nk,cells_k,np,DI,DJ,DK,...
        TI,TJ,TK,transmi);
    
    end

    fb_output = strcat(outfolder,'\fbased_',num2str(ni),'_',...
        num2str(nj),'_',num2str(nk),'.dat');
    
    if refino == 0
        
        imex_export_3d(arquivo,outfolder,fb_output,ni,nj,nk,nDI,nDJ,nDK,nPHI,...
            fbasedk_i,fbasedk_j,fbasedk_k,nPHI_include);
        
    else
        
        imex_export_3d_ref(arquivo,outfolder,fb_output,i,j,k,ni,nj,nk,cells_i,cells_j,cells_k,...
            nDI,nDJ,nDK,nPHI,PHI,KI,KJ,KK,wells,...
            fbasedk_i,fbasedk_j,fbasedk_k,nPHI_include)
        
    end
    
    %     system([imex aux1 outfolder aux2 aux fb_output aux '&'])
end